{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Construction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we build LDA, QDA, and Naive Bayes classifiers. We will demo these classes on the {doc}`wine </content/appendix/data>` dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn import datasets\n",
    "\n",
    "wine = datasets.load_wine()\n",
    "X, y = wine.data, wine.target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LDA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An implementation of linear discriminant analysis (LDA) is given below. The main method is `.fit()`. This method makes three important estimates. For each $k$, we estimate $\\pi_k$, the class prior probability. For each class we also estimate the mean of the data in that class, $\\bmu_k$. Finally, we estimate the overall covariance matrix across classes, $\\bSigma$. The formulas for these estimates are detailed in the {doc}`concept section </content/c4/concept>`.\n",
    "\n",
    "The second two methods, `.mvn_density()` and `.classify()` are for classifying new observations. `.mvn_density()` just calculates the density (up to a multiplicative constant) of a Multivariate Normal sample provided the mean vector and covariance matrix. `.classify()` actually makes the classifications for each test observation. It calculates the density for each class, $p(\\bx_n|Y_n = k)$, and multiplies this by the prior class probability, $p(Y_n = k) = \\pi_k$, to get a posterior class probability, $p(Y_n = k|\\bx_n)$. It then predicts the class with the highest posterior probability. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LDA:\n",
    "    \n",
    "    ## Fitting the model \n",
    "    def fit(self, X, y):\n",
    "        \n",
    "        ## Record info\n",
    "        self.N, self.D = X.shape\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        \n",
    "        ## Get prior probabilities \n",
    "        self.unique_y, unique_y_counts = np.unique(self.y, return_counts = True) # returns unique y and counts\n",
    "        self.pi_ks = unique_y_counts/self.N\n",
    "        \n",
    "        ## Get mu for each class and overall Sigma\n",
    "        self.mu_ks = []\n",
    "        self.Sigma = np.zeros((self.D, self.D))        \n",
    "        for i, k in enumerate(self.unique_y):\n",
    "            \n",
    "            X_k = self.X[self.y == k]\n",
    "            mu_k = X_k.mean(0).reshape(self.D, 1)\n",
    "            self.mu_ks.append(mu_k)\n",
    "\n",
    "            for x_n in X_k:\n",
    "                x_n = x_n.reshape(-1,1)\n",
    "                x_n_minus_mu_k = (x_n - mu_k)\n",
    "                self.Sigma += np.dot(x_n_minus_mu_k, x_n_minus_mu_k.T)\n",
    "            \n",
    "        self.Sigma /= self.N\n",
    "        \n",
    "        \n",
    "    ## Making classifications\n",
    "\n",
    "    def _mvn_density(self, x_n, mu_k, Sigma):\n",
    "        x_n_minus_mu_k = (x_n - mu_k)\n",
    "        density = np.exp(-(1/2)*x_n_minus_mu_k.T @ np.linalg.inv(Sigma) @ x_n_minus_mu_k)\n",
    "        return density\n",
    "            \n",
    "    def classify(self, X_test):\n",
    "        \n",
    "        y_n = np.empty(len(X_test))\n",
    "        for i, x_n in enumerate(X_test):\n",
    "            \n",
    "            x_n = x_n.reshape(-1, 1)\n",
    "            p_ks = np.empty(len(self.unique_y))\n",
    "        \n",
    "            for j, k in enumerate(self.unique_y):\n",
    "                p_x_given_y = self._mvn_density(x_n, self.mu_ks[j], self.Sigma)\n",
    "                p_y_given_x = self.pi_ks[j]*p_x_given_y\n",
    "                p_ks[j] = p_y_given_x\n",
    "            \n",
    "            y_n[i] = self.unique_y[np.argmax(p_ks)]\n",
    "        \n",
    "        return y_n\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We fit the LDA model below and classify the training observations. As the output shows, we have 100% training accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lda = LDA()\n",
    "lda.fit(X, y)\n",
    "yhat = lda.classify(X)\n",
    "np.mean(yhat == y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function below visualizes class predictions based on the input values for a model with $\\bx_n \\in \\mathbb{R}^2$. To apply this function, we build a model with only two columns from the `wine` dataset. We see that the decision boundaries are linear, as we expect from LDA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def graph_boundaries(X, model, model_title, n0 = 100, n1 = 100, figsize = (7, 5), label_every = 4):\n",
    "        \n",
    "        # Generate X for plotting \n",
    "        d0_range = np.linspace(X[:,0].min(), X[:,0].max(), n0)\n",
    "        d1_range = np.linspace(X[:,1].min(), X[:,1].max(), n1)\n",
    "        X_plot = np.array(np.meshgrid(d0_range, d1_range)).T.reshape(-1, 2)\n",
    "        \n",
    "        # Get class predictions\n",
    "        y_plot = model.classify(X_plot).astype(int)\n",
    "        \n",
    "        # Plot \n",
    "        fig, ax = plt.subplots(figsize = figsize)\n",
    "        sns.heatmap(y_plot.reshape(n0, n1).T,\n",
    "                   cmap = sns.color_palette('Pastel1', 3),\n",
    "                   cbar_kws = {'ticks':sorted(np.unique(y_plot))})\n",
    "        xticks, yticks = ax.get_xticks(), ax.get_yticks()\n",
    "        ax.set(xticks = xticks[::label_every], xticklabels = d0_range.round(2)[::label_every],\n",
    "               yticks = yticks[::label_every], yticklabels = d1_range.round(2)[::label_every])\n",
    "        ax.set(xlabel = 'X1', ylabel = 'X2', title = model_title + ' Predictions by X1 and X2')\n",
    "        ax.set_xticklabels(ax.get_xticklabels(), rotation=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAFNCAYAAACkMKB8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5wkdXnv8c9XFkVEBeUScFkxqHghsFwEPCoS0OOKVxANRI2JIJpED+ZERSVHzCEXlRjJiUGzXIQkgCFyia6wQoyIRrmsuJBFiCJBWFhYcSWKEBV4zh9VA83Qc9nZ6emamc/79erXdlX9qvr3m97pZ56nfl2VqkKSpC551LA7IEnSaAYnSVLnGJwkSZ1jcJIkdY7BSZLUOQYnSVLnGJw0qyXZIUklWdAuX5jkLVM4zqIkdyfZaPp7OeZrnpbkT2bq9YYlyU1JXjLsfmh2MTjNQ2N9WCTZL8kD7Yf03UlWJzk7yfP6tE2SG5N8Z5Kvd297zDuSfCbJZtM1nl5V9fKqOn2SfXrwZ1BVN1fVZlV1/yD6NUhJFib5cZIX9qzbvl23d7v8hiTfSHJPkkuG1tlRkpyTZOmodecn+WT7/BVJvp7kriS3JzkpyeOH01vNJIOTRrutqjYDHg/sA1wPfC3JAaPa7QtsDfxqv+DVx6va4+4OPA/4o9EN2oDn/8n1VFWrgaOBk5Ns0q7+W+AzVXV5u7wOOAH4yBC6OJ7fB16X5NcBkvwGsBvw/nb7E4E/AbYDng0sBI4fQj81w/wgUF/VWF1VHwJOBj46qslbgH8GLmifT/a4twIXAjsDJLkkyZ8m+TfgHppg98QkpyRZk+TWJH8yUm5LslGSv0hyZ5IbgVf0Hr893hE9y29Lcl2Snyb5TpLdk/w9sAj4QpvNva9PeXC7JJ9Psi7JDUne1nPMD7cZ5d+1x702yZ49249u+/3TJP/RJ7D32jLJxW3bryZ5anuMv0ny8VFj+0KSd49xnJOANcCxbVlzJ3r+AKiqf6mqs4HbxunLyOtskWRZkh+22deyJAtH/YyPS/Jvbb8vSrJlz/Y3J/lBkh8lOWa816qq24E/BE5Ksgj4f8Dbq+rudvuZVbW8qu6pqh+343zBRGPQHFBVPubZA7gJeEmf9fsBq/us3x94AHhcu7wp8BPgQOB1wJ3AoyfzesD2wLXAce3yJcDNwHOBBcDGwPk0f/k/jiY7u4LmAwvgHTTZ3PbAk4CvAAUs6DneEe3z1wO30mRqAZ4OPLXfzwDYYdRxvgqcCGwCLAZ+CBzQbvsw8N/t+DcC/hy4rN22E3ALsF3PcXcc4+dyGvBTmiz0McBfAV9vt+1FE0ge1S5vSRO8txnn57wj8F/Aj4H9x2hzBHDJBP8/nty+r5vSZND/BJzfs/0S4PvAM4HHtssfabc9B7i7Z0x/CdxHn/9vo17zS+3/o9MnaHcC8Nlh/w75GPzDzEmTcRvNh/vm7fLBwM+Bi4BlNEHlFf13fdD5Se4Cvk7zwf9nPdtOq6prq+o+moDzcuDdVfWzqloLfAI4tG37BuCEqrqlqtbRBIaxHAF8rKqurMYNVfWDiQabZHvghcDRVfXfVbWSJnt8c0+zr1fVBdWco/p7YNd2/f00H8rPSbJxVd1UVd8f5+W+WFWXVtXPgWOA5yfZvqquoAk0I1nXoTRB5Y5xjvUDmvfqJ8ClE41zLFX1o6o6p5ps5afAnwIvHtXsM1X13aq6FzibJoADHAIs6xnT/6H5w2YiX6MJiv8wVoMkL6XJ0j+0fiPSbGRw0mQ8hSaruKtdfgtwdlXd134AncvEpb3XVtXmVfXUqvq99kNtxC09z59Kkz2taU+C30WTRW3dbt9uVPvxgs32NH/hr6/tgHXtB3Pv6zylZ/n2nuf3AJskWVBVNwDvpsmu1ib5bJLtxnmtB8dSTSlrXfv6AKcDb2qfv4kmCI7n/cCPgLXAeyZoO6Ykmyb527Y0NxLoNs/DZzKOHv/IBJeHvT9V9bO2T+O93jPa/p4IfDzJxn3a7AOcCRxSVd+dwrA0yxicNBkHAVdV1c/acw/7A29qZ0/dTvPX8oG95x3WU++l8W+hycq2bIPZ5lX1hKp6brt9DU3QGbFonOPeQlPqmug1R7sNeNKoWWGLaEqEE6rmPMkLaQJt8cjzdb0eHEs7g/FJPHRe6B+A1yTZlWYywPljHSTJc4D30mSLhwMfbD/0p+IPacqTe1fVE2hKdNBkzxN52PuTZFOajGisfocmKz0BeBfwM5rJHb1tdgM+D7y1qr48+WFoNjM4zV8bJ9mk57Ggd2M7c+4pSY6l+cD7YLvpzcB3aT68FrePZwKrgcM2tFNVtYamXPjxJE9I8qgkOyYZKSudDfyvNNOnt+ChWV39nAy8J8ke7XiePjLhALgD+NUx+nAL8A3gz9ufzS40H/hnTNT/JDsl2T/JY2jOS91LU+oby4FJXpjk0cBxwOXt61PNLLwraTKmc0Zlm72v+SjgFJoS5vVVdQ3NxIKl7Yf/yESSTWhKsI9qx/WIDKX1+LbfdyV5EnDsROPu8TnglT1j+r+M/znzuzTn0/6sqh6g+Tm/L8mz2n7vDCwH3lVVX1iPfmiWMzjNXxfQfACNPD7crt8uyd00J7WvBH4N2K+qLmq3vwU4sapu730An2Y9Zu1N4LeARwPfoTm5/zlg23bbSTQnz68GrqIpKfZVVf9Ec77kTJqJB+fTZCbQnKv6o7Z02K8EdhjNZIbbgPOAY6vq4kn0/TE007XvpCl9bc1Dgb2fM2k+/NcBewBvHLX9dJr3YLyS3lE0kxc+1rPuOOBXaP6wgOaPinuBTwEvap+fNMbxTqCZ6HAncBlNcJiUqrqWZnr4mTRZ1I9p/nB5hPbc3p8Bh1fVL9r9vwN8nGb2XmiyuK2AU/LQ9++unWx/NHulypsNSl2VZF+a8t4ObWYhzQtmTlJHtWW3o4CTDUyabwxOUgcleTbN7Mhtacps0qyV5nJaX0nzhfhrkxw14T6W9SRJg5RkW2DbqrqqnQX7LZqvl4x5bc6hZE5JTk2yNsmqnnVPai/j8r323y3G2HdRe7mU69JcjmaHmeq3JGn9VdWaqrqqff5T4Doe/r3BRxhWWe80YMmode8HvlxVzwC+zNhThP8OOL6qnk1ziZe1g+qkJGl6tQnFbsDl47VbMN7GQamqS/tkPK+hubYbNNNnL+GRX8Z7Ds21zy5uj3P3ZF5v2VWrO127XLLmymF3QdIsteAVB03my9GT9q07vzGlz8s9t3rB24Eje1YtrarRt0PZDDiH5vJkPxnveF2aELFN+wXMkS9ibt2nzTNpvhh4bpJvJzk+Y9wcLsmRSVYkWbH83Am/OylJ2gBVtbSq9ux5jA5MG9MEpjOqaszvJ44YSua0ARbQfIFwN5orWf8j8Ns0345/mPYHsxS6nzkt3/ah2yGZRUmaa9ovVJ8CXFdVfzmZfbqUOd3RzugYmdnR71zSauDbVXVjewXr82luXidJ6q4X0FylZP8kK9vHgePt0KXM6fM0l7/5CA/dyG60K4EtkmxVVT+kuQDpipnroiRpfVXV15nchYMfNJTglOQsmskPWyZZTXNtsY8AZyc5nKZk9/q27Z7AO6rqiKq6v70O2pfbNPFbjH19sFnJEp8kDW+23lhXr37E7ayragUPXbySdqbeLgPqmiSpA7p0zkmSJKBb55w0iiU+SfOVmZMkqXPMnGaJ3iwKzKQkzW1mTpKkzjE4SZI6x7LeLOVkCUlzmZmTJKlzDE6SpM6xrDcHWOKTNNeYOUmSOsfgJEnqHMt6c4wlPklzgZmTJKlzzJzmMLMoSbOVmZMkqXMMTpKkzrGsN09Y4pM0m5g5SZI6x+AkSeocy3rzkCU+SV1n5iRJ6hwzp3nOLEpSF5k5SZI6x+AkSeocy3p6UG+JDyzzSRoeMydJUucYnCRJnWNZT2NyJp+kYTFzkiR1jsFJktQ5lvU0KZb4JM0kMydJUueYOWm9mUVJGrSBZk5JTk2yNsmqnnWvT3JtkgeS7DnGfpskuSLJ1W3bP+7ZdkCSq5KsTPL1JE8f5BgkSTNv0GW904Alo9atAg4GLh1nv58D+1fVrsBiYEmSfdptnwLeWFWLgTOBP5rWHkuShm6gZb2qujTJDqPWXQeQZLz9Cri7Xdy4fdTIZuAJ7fMnArdNW4e13izxSRqEzp5zSrIR8C3g6cDfVNXl7aYjgAuS3Av8BNhnjENIkmapzs7Wq6r729LdQmCvJDu3m/4AOLCqFgKfAf6y3/5JjkyyIsmK5eeeMTOdliRNi85mTiOq6q4kl9Ccd7oD2LUni/pHYPkY+y0FlgIsu2p19Wuj6WWJT9J06WTmlGSrJJu3zx8LvAS4Hvgx8MQkz2ybvhS4bji9lCQNykAzpyRnAfsBWyZZDRwLrAP+GtgK+GKSlVX1siTbASdX1YHAtsDp7XmnRwFnV9Wy9phvA85J8gBNsHrrIMegqTGLkrQhBj1b77AxNp3Xp+1twIHt82uA3cY45nn99pckzR2dLOtJkua3zk+I0OxniU/S+jJzkiR1jsFJktQ5lvU0o3pLfGCZT1J/Zk6SpM4xOEmSOseynobKmXyS+jFzkiR1jpmTOsMsStIIMydJUucYnCRJnWNZT51kiU+a38ycJEmdY3CSJHWOZT11niU+af4xc5IkdY6Zk2YVsyhpfjBzkiR1jsFJktQ5lvU0a1nik+YuMydJUucYnCRJnWNZT3OCJT5pbjFzkiR1jsFJktQ5lvU05/SW+MAynzQbmTlJkjrHzElznpMlpNnHzEmS1DkGJ0lS51jW07xiiU+aHcycJEmdY3CSJHWOZT3NW5b4pO4aWOaU5NQka5Os6ll3fJLrk1yT5Lwkm09233b9PyZZ2T5uSrJyUP2XJA3PIMt6pwFLRq27GNi5qnYBvgt8YD32pap+o6oWV9Vi4Bzg3Gnrrea15ds+78GHpOEbWHCqqkuBdaPWXVRV97WLlwELJ7tvryQB3gCcNT29lSR1yTAnRLwVuHCK+74IuKOqvjdWgyRHJlmRZMXyc8+Y4stIkoZhKBMikhwD3AdMNWocxgRZU1UtBZYCLLtqdU3xdTQPOVFCGr4ZD05J3gK8EjigqtY7aCRZABwM7DHdfZMkdcOMBqckS4CjgRdX1T1TPMxLgOuravX09UySBLDm5kVT23HL6e3HIKeSnwV8E9gpyeokhwOfBB4PXNxOB/9023a7JBdMsO+IQ3EihGaIs/ik4RhY5lRVh/VZfcoYbW8DDpxg35Ftv73BnZMkdZqXL5IkdY6XL5ImyVl80swxc5IkdY6ZkzQFoydImElJ08vMSZLUOQYnSVLnWNaTpoGTJaTpZeYkSeocg5MkqXMs60nTzBKftOHMnCRJnWNwkiR1jmU9aYAs8UlTY+YkSeocMydphphFSZNn5iRJ6hyDkySpcyzrSUNgiU8an5mTJKlzDE6SpM6xrCcNmSU+6ZHMnCRJnWPmJHWIWZTUMHOSJHWOwUmS1DmW9aSO6i3xgWU+zS9mTpKkzjE4SZI6x7KeNEtcvfc2Dz7f9fI7htgTafDMnCRJnWNwkiR1jmU9aZZYc/Oihxb2fuipJT7NRWZOkqTOMXOSZrneiRK9zKg0mw0sc0pyapK1SVb1rDs+yfVJrklyXpLNx9j3qCSrklyb5N096z+c5NYkK9vHgYPqvyRpeAZZ1jsNWDJq3cXAzlW1C/Bd4AOjd0qyM/A2YC9gV+CVSZ7R0+QTVbW4fVwwkJ5LkoZqYMGpqi4F1o1ad1FV3dcuXgYs7LPrs4HLquqetu1XgYMG1U9prrp6720efEizzTAnRLwVuLDP+lXAvkmenGRT4EBg+57t72zLgqcm2WImOipJmlnjBqckT0iyY5/1u2zIiyY5BrgPOGP0tqq6DvgoTQlwOXB12xbgU8COwGJgDfDxcV7jyCQrkqxYfu4jXkaSNIP6zUMYz5jBKckbgOuBc9qJCb2XSD5tAzr4FuCVwBurqvq1qapTqmr3qtqXpjT4vXb9HVV1f1U9AJxEc16qr6paWlV7VtWeSw5+41S7K80JlvjUAafxyHkIYxovc/ogsEdVLQZ+B/j7JAe32zKVniVZAhwNvLqq7hmn3dbtv4uAg4Gz2uVte5odRFMClCR1XL95COMZ73tOC6pqTXvQK5L8OrAsyUKgb8bTK8lZwH7AlklWA8fSzM57DHBxEmgmPrwjyXbAyVU1MjX8nCRPBn4J/H5V/bhd/7Eki9vXvwl4+2QHKqkxXvbkd6PUFeMFp58k2bGqvg9QVWuS7AecDzx3ogNX1WF9Vp8yRtvbaCY+jCy/aIx2b57odSVJMy/JkcCRPauWVtXSqR5vvOB0NKPKd1X107Y094jvJ0mS5q82EE05GI023jmn04HXJXkwgCXZBvgM8Krp6oCk7nDihLpivOC0B/A04NtJ9k9yFHAF8E0edk1kSZLG185D+CawU5LVSQ4fr/2YZb12EsI72qD0L8BtwD5VtXo6OyxJmvvGmIcwpjGDU3tR1o/SZElLaCYsXJjkqKr61w3qpaQN0ntvp20X3TyQ1/C28Bqm8SZEXAWcSDOV+z7gonYa94lJfrC+UVCSpMkaLzjtO7qEV1Urgf+R5G2D7ZYkaT4b75zTmOeWquqkwXRHUhd5Q0PNNG/TLknqHG/TLmnKnDShQTFzkiR1jsFJktQ5lvUkTQsnTWg6mTlJkjrH4CRJ6hzLepIGanS5zzKfJsPMSZLUOWZOkmaU343SZJg5SZI6x+AkSeocy3rSLDcT93YaFL8bpbGYOUmSOsfgJEnqHMt6kjrHGX0yc5IkdY7BSZLUOZb1JHWaJb75ycxJktQ5Zk6SZg2/FzV/mDlJkjrH4CRJ6hzLepJmPe8ZNfeYOUmSOsfgJEnqHMt6kuYcZ/XNfgPNnJKcmmRtklU9645Lck2SlUkuSrJdn/1+vd0+8vjvJK8d1eavk9w9yP5LkoZj0GW904Alo9YdX1W7VNViYBnwodE7VdVXqmpx22Z/4B7gopHtSfYENh9Yr6VZas3Nix586JGu3nubBx/qtoEGp6q6FFg3at1PehYfB9QEhzkEuLCq7gFIshFwPPC+aeyqJKlDhjIhIsmfJrkFeCN9MqdRDgXO6ll+J/D5qlozwWscmWRFkhXLzz1jwzosSZpRQ5kQUVXHAMck+QBNsDm2X7sk2wK/BnypXd4OeD2w3yReYymwFGDZVasnys4kzTNeULbbhj2V/EzgdeNsfwNwXlX9sl3eDXg6cEOSm4BNk9ww2C5KkmbajAenJM/oWXw1cP04zQ+jp6RXVV+sql+pqh2qagfgnqp6+mB6KkkaloGW9ZKcRVOC2zLJapry3YFJdgIeAH4AvKNtuyfwjqo6ol3eAdge+Oog+yhJfi+qewYanKrqsD6rTxmj7QrgiJ7lm4CnTHD8zTakf5Kkbhr2OSdJkh7ByxdJ0hjm44y+JWuunOKeC6e1H2ZOkqTOMXOSpEnwnlEzy8xJktQ5BidJUudY1pOkKfC7UYNl5iRJ6hyDkySpcyzrSXNU7w0Ht1108xB7Mr/Mx+9GDYKZkySpc8ycJGlAzKKmzsxJktQ5BidJUudY1pOkGeD3otaPmZMkqXMMTpKkzrGsJ0lD5Iy+/sycJEmdY3CSJHWOZT1J6oixZvTB/Cv5mTlJkjrHzEmSZoH5NnHCzEmS1DkGJ0lS51jWk+YB7+00t/SW+PYYYj8GycxJktQ5BidJUucYnCRJnWNwkiR1jsFJktQ5BidJUucYnCRJnWNwkiR1zsCCU5JTk6xNsqrPtvckqSRbjrHv/UlWto/P96w/I8l/JFnVHn/jQfVfkjQ8g8ycTgOWjF6ZZHvgpcB4X1O/t6oWt49X96w/A3gW8GvAY4Ejpq+7kqSuGFhwqqpLgXV9Nn0CeB9QUzjmBdUCrgAWblgvJUldNKPnnJK8Gri1qq6eoOkmSVYkuSzJa/scZ2PgzcDyQfRTkjRcMxackmwKHAN8aBLNF1XVnsBvAick2XHU9hOBS6vqa+O83pFtgFux/NwzptxvSdLMm8mrku8IPA24Ogk0JbmrkuxVVbf3Nqyq29p/b0xyCbAb8H2AJMcCWwFvH+/FqmopsBRg2VWr17uEKEkanhkLTlX178DWI8tJbgL2rKo7e9sl2QK4p6p+3s7mewHwsXbbEcDLgAOq6oGZ6rskaWYNcir5WcA3gZ2SrE5y+Dht90xycrv4bGBFkquBrwAfqarvtNs+DWwDfLOdZj6ZEqEkaZYZWOZUVYdNsH2HnucraKeFV9U3aKaK99vHmyNK0jzgFSIkSZ1jJiLNM96yXbOBmZMkqXMMTpKkzjE4SZI6x+AkSeocg5MkqXMMTpKkzjE4SZI6x+AkSeocg5MkqXMMTpKkzjE4SZI6x+AkSeocg5MkqXMMTpKkzjE4SZI6x/s5SfNY772dwPs7qTvMnCRJnWNwkiR1jsFJktQ5BidJUucYnCRJnWNwkiR1jsFJktQ5BidJUucYnCRJnWNwkiR1jsFJktQ5BidJUucYnCRJnWNwkiR1jsFJktQ5BidJUucMJTglOTXJ2iSretZ9OMmtSVa2jwP77LdJkiuSXJ3k2iR/PLM9l+a2NTcvevAhDdOwMqfTgCV91n+iqha3jwv6bP85sH9V7QosBpYk2WeA/ZQkDcFQglNVXQqsm8J+VVV3t4sbt4+azr5Jkoava+ec3pnkmrbst0W/Bkk2SrISWAtcXFWXz2wXJUmD1qXg9ClgR5py3Rrg4/0aVdX9VbUYWAjslWTnfu2SHJlkRZIVy889Y1B9liQNwIJhd2BEVd0x8jzJScCyCdrfleQSmnNXq/psXwosBVh21WpLf5I0i3Qmc0qybc/iQfQJOEm2SrJ5+/yxwEuA62emh5KkmTKUzCnJWcB+wJZJVgPHAvslWUwzweEm4O1t2+2Ak6vqQGBb4PQkG9EE1rOratwMS5I0+wwlOFXVYX1WnzJG29uAA9vn1wC7DbBrkqQO6ExZT5KkEQYnSdLAJVmS5D+S3JDk/RO1NzhJkgaqnSfwN8DLgecAhyV5znj7GJwkSYO2F3BDVd1YVb8APgu8ZrwdDE6SpEF7CnBLz/Lqdt2YOvMl3EF65e4LA81VI9ov585qjqNb5sI4+o9h4XA6swHmwnsBwx3HglcclKnsl+RI4MieVUt7xtDvmONeHGG+ZU5HTtxkVnAc3TIXxjEXxgCOY2iqamlV7dnz6A2uq4Hte5YXAreNd7z5FpwkSTPvSuAZSZ6W5NHAocDnx9thXpT1JEnDU1X3JXkn8CVgI+DUqrp2vH3mW3Ca9bXoluPolrkwjrkwBnAcndXeQLbfTWT7SpUX7JYkdYvnnCRJnTMnglN759y1SR5xm412+2vaO+yubG9A+MKebYuSXJTkuiTfSbLDTPW7Tz/HHUdPu+cluT/JIaPWPyHJrUk+Odiejtu3id6L/ZL8V/terEzyocnuO5M2cBxHJVmV5Nok7565Xvft54Q/03YsK9v+frVn/eZJPpfk+vb34/kz0+u+fZzo/Xhvz3uxqv39eFLP9o2SfDvJUO9iMIlxPDHJF5Jc3b4fv9OzbXmSu4Y9hhlTVbP+AewL7A6sGmP7ZjxUwtwFuL5n2yXAS3vabdrVcbRtNgL+laZ2e8iobX8FnAl8sqtjoLlVyrKpjr/r4wB2prkX2aY053T/BXhGh8exOfAdYFG7vHXPttOBI9rnjwY27+o4RrV9FfCvo9b97/Z3o+//va6MA/gg8NH2+VbAOuDR7fIB7diGOoaZesyJzKmqLqV5E8fafne17y7wONovf7XXdlpQVRf3tLtn0P0dy0TjaL0LOAdY27syyR7ANsBFg+nd5ExyDNO+73TbgL48G7isqu6pqvuAr9LcPHMoJjGO3wTOraqb2/ZrocnCaT5IT2nX/6Kq7hpwd8e0nu/HYcBZIwtJFgKvAE4eQNfWyyTGUcDjk4Tmj+V1wH3tvl8GfjrwTnbEnAhOk5HkoCTXA18E3tqufiZwV5Jz25T/+PYChZ2U5Ck0H3SfHrX+UcDHgfcOo19T8Py2bHFhkucOuzMboN84VgH7Jnlykk1p7kW2/diHGLpnAlskuSTJt5L8Vrv+V4EfAp9pfzdOTvK44XVzctqf+RKaP+BGnAC8D3hgKJ1aP5+k+QPnNuDfgaOqajb0e9rNm+BUVedV1bOA1wLHtasXAC8C3gM8j+YX8reH0sHJOQE4uqruH7X+94ALquqWPvt0zVXAU6tqV+CvgfOH3J+p6juOqroO+ChwMbAcuJr2L9+OWgDsQZNZvAz4P0me2a7fHfhUVe0G/AyY8DYHHfAq4N+qah1AklcCa6vqW8Pt1qS9DFgJbAcsBj7ZZrHzzrwJTiPatHrHJFvSXFLj29VcKfc+mg+Y3YfawfHtCXw2yU3AIcCJSV4LPB94Z7v+L4DfSvKRofVyHFX1k6q6u31+AbBx+17MKuONo6pOqardq2pfmrLM94bY1YmsBpZX1c+q6k7gUmDXdv3qqrq8bfc5uv27MeJQekp6wAuAV7e/G58F9k/yD8Po2CT9Dk2ZtarqBuA/gWcNuU9DMS+CU5KntzVckuxOc3L3RzSX1NgiyVZt0/1pTg53UlU9rap2qKodaD4sfq+qzq+qN1bVonb9e4C/q6pO/pWb5Fd63ou9aP4P/mi4vVp/440jydbtv4uAg3n4h2XX/DPwoiQL2pLY3sB1VXU7cEuSndp2B9Dh3w1oZroBL6YZEwBV9YGqWtj+bhxKM1HiTUPq4mTcTPOzJsk2wE7AjUPt0ZDMiStEJDmLZvbUlklWA8cCGwNU1aeB19FkE78E7gV+o50gcX+S9wBfbj9ovgWcNIQhAJMaR+dNYgyHAL+b5D6a9+LQkckq/fatqlNmfBBj9IVJjgM4J8mTgV8Cv19VP57p/o+YaBxVdV2S5cA1NOdkTq6qkWnO7wLOSHMttBtp/qofikn+bhwEXFRVPxtKJydhEuM4Djgtyb/TXMn76DajJcnXaLKozdp9D6+qL838KGaGV4iQJHXOvCjrSZJmF4OTJKlzDE6SpM4xOEmSOsfgJEnqHIOTNEVJtk/ynyNXv06yRbv81Hl3BWlpmhmcpClqL2ICEs8AAADWSURBVBf1KWDkahwfAZZW1Q+A44E3D6tv0mxncJI2zCeAfdLct+mFNBfgnXdXkJam25y4QoQ0LFX1yyTvpbnI6/+sql8Mu0/SXGDmJG24lwNraG40KGkaGJykDZBkMfBSYB/gD5JsO+QuSXOCwUmaovZiwZ8C3t3eSfZ4mluWSNpABidp6t4G3FxVF7fLJwLPSvLi9grS/wQckGR1kpcNrZfSLORVySVJnWPmJEnqHIOTJKlzDE6SpM4xOEmSOsfgJEnqHIOTJKlzDE6SpM4xOEmSOuf/AxSEP00yNPpcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X_2d = X.copy()[:,2:4]\n",
    "lda_2d = LDA()\n",
    "lda_2d.fit(X_2d, y)\n",
    "graph_boundaries(X_2d, lda_2d, 'LDA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QDA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The QDA model is implemented below. It is nearly identical to LDA except the covariance matrices $\\bSigma_k$ are estimated separately. Again see the {doc}`concept section </content/c4/concept>` for details. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "class QDA:\n",
    "    \n",
    "    ## Fitting the model\n",
    "    \n",
    "    def fit(self, X, y):\n",
    "        \n",
    "        ## Record info\n",
    "        self.N, self.D = X.shape\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        \n",
    "        \n",
    "        ## Get prior probabilities \n",
    "        self.unique_y, unique_y_counts = np.unique(self.y, return_counts = True) # returns unique y and counts\n",
    "        self.pi_ks = unique_y_counts/self.N\n",
    "        \n",
    "        \n",
    "        ## Get mu and Sigma for each class\n",
    "        self.mu_ks = []\n",
    "        self.Sigma_ks = []\n",
    "        for i, k in enumerate(self.unique_y):\n",
    "            \n",
    "            X_k = self.X[self.y == k]\n",
    "            mu_k = X_k.mean(0).reshape(self.D, 1)\n",
    "            self.mu_ks.append(mu_k)\n",
    "            \n",
    "            Sigma_k = np.zeros((self.D, self.D))\n",
    "            for x_n in X_k:\n",
    "                x_n = x_n.reshape(-1,1)\n",
    "                x_n_minus_mu_k = (x_n - mu_k)\n",
    "                Sigma_k += np.dot(x_n_minus_mu_k, x_n_minus_mu_k.T)\n",
    "            self.Sigma_ks.append(Sigma_k/len(X_k))\n",
    "     \n",
    "    ## Making classifications \n",
    "    \n",
    "    def _mvn_density(self, x_n, mu_k, Sigma_k):\n",
    "        x_n_minus_mu_k = (x_n - mu_k)\n",
    "        density = np.linalg.det(Sigma_k)**(-1/2) * np.exp(-(1/2)*x_n_minus_mu_k.T @ np.linalg.inv(Sigma_k) @ x_n_minus_mu_k)\n",
    "        return density\n",
    "    \n",
    "    def classify(self, X_test):\n",
    "        \n",
    "        y_n = np.empty(len(X_test))\n",
    "        for i, x_n in enumerate(X_test):\n",
    "            \n",
    "            x_n = x_n.reshape(-1, 1)\n",
    "            p_ks = np.empty(len(self.unique_y))\n",
    "        \n",
    "            for j, k in enumerate(self.unique_y):\n",
    "\n",
    "                p_x_given_y = self._mvn_density(x_n, self.mu_ks[j], self.Sigma_ks[j])\n",
    "                p_y_given_x = self.pi_ks[j]*p_x_given_y\n",
    "                p_ks[j] = p_y_given_x\n",
    "            \n",
    "            y_n[i] = self.unique_y[np.argmax(p_ks)]\n",
    "        \n",
    "        return y_n\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9943820224719101"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qda = QDA()\n",
    "qda.fit(X, y)\n",
    "yhat = qda.classify(X)\n",
    "np.mean(yhat == y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below plot shows predictions based on the input variables for the QDA model. As expected, the decision boundaries are quadratic, rather than linear. We also see that the area corresponding to class 2 is much smaller than the other areas. This suggests that either there were fewer observations in class 2 or the estimated variance of the input variables for observations in class 2 was smaller than the variance for observations in other classes ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAFNCAYAAACkMKB8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5gkZXn38e9PFiOICgoSjq4BQQ3qioBnJKJxRTygREE8i6jRBPPGs4mivnmjEg+JCGY5CCqiJuIhCBuIisQEDwsuuAiKUYSVVURARIxh4X7/qBrSDj09s7PT3TUz38919bVdVU9V30/3Tt99P/V0daoKSZK65E7jDkCSpMlMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzjE5ad5L8qIkX+tZvinJH8ziOIclOXtuo5v2Ma9I8oRRPuaoJVmapJIsGXcsmj9MTotM+0b+nSQ3J/lpkmOT3KNn+1FJbknyq/b2/STHJNmuz7Hum+S2JMdO85gTb043tbcrkrxxGP0DqKotquqHM4xpSc9+p1bVHw8rrmFK8tT29bxnz7qnJ/nJxOub5J3ta78+yVFjC7ZHki3a/w/P7Vl3tyRXJjm4XX5dkjXt/8cfJXnd+CLWqJicFpEkfwm8G3gdcA/gEcBS4Owkm/Y0/VRV3Q24J3AQ8PvABX0S1AuA64FDkvzeDELYsqq2AA4F3ppkeZ8Y/XQ9C1X1L8CXgfcDJNkSOA54ZVX9sm32A+D1wBfHEmQfVXUTcATw90m2aVe/B1hVVf/cLofm/9pWwHLg1UkOGXmwGimT0yKR5O7A24E/q6qVVXVLVV0BPBu4L/Dcyfu0bS4BngP8HPjLSU1eAPwVcAvw1JnGUlXnA5cAe7SxVZJXJbkcuLxdd/8k5yS5Lsn3kjy7py/3SvKFJDcm+Sawy6S+VpJd2/ubJXlvkh8n+WWSryXZDDivbX5DW809ss/w4KOSfKvd71tJHtWz7dy2EvmP9hP92Um2brfdJcnHk/wiyQ3tvtsOeEr2TvLdJNcn+UiSu7THWZPk9uc1yaZJrk2ybIrj/Dnw5CRPoklSX62qL/Q876dU1VnArwbEMvFY+yQ5v41/XVs933nSc/yKJJe3cX8oSdptmyT5uzbWHwJPGfRYVXU2TcL8hyT70fyffFXP9vdU1YVVtb6qvgd8Hnj0dH3QPFdV3hbBjeYT53pgSZ9tpwCntvePAj7ep807gG/0LD8W+C3Np9kPAl8Y8NhLgQKW0HwKfjRwM7B/u72Ac2gqtc2AuwJXAS9u99kTuBb4w7b9J4FPt+32AH4CfK3n8QrYtb3/IeBcYAdgE+BRwO/1xtSz34smjtPGcj3w/DaGQ9vle7XbzwX+C9itjflc4F3ttpcD/wJs3j7mw4C7T/HcXAGsAXZqH/M/gP/bbns9TRU70fbpwHemeZ0PbZ+rnwPbTNHm48BR0xznYTSV9ZL2uboUeM2k5/gMYEtg5/bxlrfbXgFc1tOnr0x+rvs83lbAujb2Fw9oF+DbwCvG/Tflbbg3K6fFY2vg2qpa32fbOmCbPut7XU3zRjPhhcBZVXU98AmaT+z3nuYY1wLXAScAb6yqL/Vs+9uquq6qfgMcCFxRVR+p5tPyhcBngIOTbAI8C3hrVf26qtbQJNc7SHIn4CXAkVX1k6q6tar+s6p+O02c0Hzav7yqPtbGcBrNG25vhfiRqvp+G/OngYmK5hbgXjQJ8taquqCqbhzwWMdU1VVVdR3wNzQJBpokckBb9UKTKD82TdxfpxmyPbuqfj6DfvbVxvz1tu9XAP8IPG5Ss3dV1Q1VdSVNApro/7OBD/T06W9n8HjX01TTmwOnD2h6FM2Iz0c2pD+af0xOi8e1wNZTnNPZjuaT7yA70CQW2mGxPwFOhduH6a6kz9DgJFtX1VZV9YCq+odJ267quX8f4OHtkNINSW4ADqM597UNzaf53vY/nurxgLvQVDgbavs+x/0xzfMw4ac9928Gtmjvfwz4V+CTSa5O8p5J5/Qmm9yX7QGq6mqaSupZ7TmkJ9M+5wOsAD5Kk9QeNU3bKSXZLckZ7SSLG4H/R/N89pqq/9szs9en9/GeR1Oh/RvNedF+bV5NM5T8lBl+wNA8ZnJaPM6nGYZ7Zu/KJHeledP76lQ7thXIU4F/b1cdBNwdOLZ98/opzZv2CzYivt7L419Fc75ky57bFlX1Spokup5myGjCzlMc81rgv5l0TqrP4/VzNU2S7LUzzRDiQNWcq3t7VT2QZhjxQAY/N5P7cnXP8inA82g+DJxfVVM+fpKXtsf6U+DNwPG954k20HE0leL9quru7fEyw33XMbPXB4C24n4/8DKaIdFnJ9l3UpuXAG+kGQpeO8M4NI+ZnBaJamZsvR34YJLl7cn1pcA/0byJ3+ETedvmAcBpNFXL+9pNLwROAh5EM5SzjOY80rIkD5qDcM8Adkvy/DaGTZPsneQBVXUrzbDPUUk2T/LANp5+fb6tjfN9SbZvT9Q/sp1Z+HPgNmCq70Od2cbw3CRLkjwHeGAb20BJ/ijJg9ohyBtphvluHbDLq5LsmGYa+JuBT/Vs+xzNObcjaSqiqR5ze+Bo4GVtVfFh4BfAW3rabNpOtrgTsKSduLHJFIe8Wxv7TUnuD7xyYKd/16eBP2/7tBVNUhnkGOBzVfWVqlpHc67t+PZ1IslhNJXbE2uarwhoARn3SS9vo70BL6U5Af/fNNXDucD2PduPonkzvQn4Nc3suWOBHdrtO9BULg/qc+wzgb/rs34pA06I0zOBoWfd7jQzuH5O8yb7ZWBZu20bmiRxI/BN4J1MPSFiM+ADNBXPL2lm6W3WbntHe/wbaE7+v2jScR4DXNDudwHwmJ5t5wKH9yzfvi/NOaPvtc/fz4B/GND3K4A3Ad9t4zgF2HxSmxPaY20x4HX9HHBsn+fwl/zvRJKT2+em9/aiKY63L03ldBNNxfyOqZ7jnmNPTORYQlMJ/QL4Ec3Mu76vP/AMmkpxy0nrvwT8TXv/R/zv/8mJ24fH/bfkbbi3tC++FqF2qOTtwKOrOamtDkryVmC3qnreuGORRsUvPC5iVXVSkltozouYnDqoHep7Kc1MPWnR8JzTIlfNVOlPjjsO3VGSl9FMDjmrqs6brr3UVUl2SvKVJJcmuSTJkdPu47CeJGmY0lz6bLuqujDJ3WjO4T6jqr471T5jqZySnJTkmiRretbdM83lai5v/91qin13bi8Vc2l7yZelo4pbkrThqmpdNV+mp6p+RXPFkR0G7TOuYb2TaS6n0+uNwJeq6n40M3Wmmn76UeDoqnoAsA9wzbCClCTNrbageCjwjUHtxjIhoqrO61PxPB3Yr71/Cs1U3Tf0Nmi/07Kkqs5pj3PTTB7vjAvXOnY5A8vXfWvcIUjaQEuectBMvxw9Ixdc+5+zer/ca5tHv5zmCvMTVlTVit42SbaguRTZa2rwJb06NSFi22q+gEf7b7/rtO1GcxXp05N8O8nRU32JMMkRSVYlWbXy9Omu+CJJ2hhVtaKq9uq5TU5Mm9IkplOratD1E4H5N5V8Cc3VsB9KM/X5UzRffjxxcsP2iVkBVk4ztXK7vX9n2UpK0lxof07lRODSqnrfdO2hW5XTz9oZHRMzO/qdS1oLfLuqfljN1bUnLu0iSequR9N8V+/xSVa3twMG7dClyukLNNdIe1f77+f7tPkWsFWSbar5OYDHA6tGF6IkaUNV1deY+YWDgfFNJT+N5irZuydZ215N+V3AE9P8GuoT22WS7JXkBIBqLvr5WuBLSb5D09njx9EHSdLwjGu23qFTbNq/T9tVwOE9y+cADx5SaJKkDujSOSdJkgCTkySpg0xOkqTO6dJsPXVM7/ee/M6TpFGycpIkdY7JSZLUOSYnSVLnmJwkSZ1jcpIkdY6z9TQjztyTNEpWTpKkzjE5SZI6x+QkSeock5MkqXOcEKEN5uQIScNm5SRJ6hyTkySpc0xOkqTOMTlJkjrH5CRJ6hxn62mjOHNP0jBYOUmSOsfKSXPGKkrSXLFykiR1jslJktQ5JidJUueYnCRJnWNykiR1jrP1NBTO3JO0MaycJEmdY3KSJHWOyUmS1DkmJ0lS5zghQkPn5AhJG2qolVOSk5Jck2RNz7o/SXJJktuS7DXFfndJ8s0kF7Vt396zbf8kFyZZneRrSXYdZh8kSaM37GG9k4Hlk9atAZ4JnDdgv98Cj6+qhwDLgOVJHtFuOw44rKqWAZ8A/mpOI5Ykjd1Qh/Wq6rwkSyetuxQgyaD9CripXdy0vdXEZuDu7f17AFfPWcCSpE7o7DmnJJsAFwC7Ah+qqm+0mw4HzkzyG+BG4BFTHEKSNE91drZeVd3aDt3tCOyTZI92018AB1TVjsBHgPf12z/JEUlWJVm18vRTRxO0JGlOdLZymlBVNyQ5l+a808+Ah/RUUZ8CVk6x3wpgBcAZF66tfm00er0z98DZe5L662TllGSbJFu29zcDngBcBlwP3CPJbm3TJwKXjidKSdKwDLVySnIasB+wdZK1wNuA64APAtsAX0yyuqqelGR74ISqOgDYDjilPe90J+DTVXVGe8yXAZ9JchtNsnrJMPug4fI7UJL6GfZsvUOn2PTZPm2vBg5o718MPHSKY3623/6SpIWjk8N6kqTFzeQkSeock5MkqXNMTpKkzun895y0eDhzT9IEKydJUueYnCRJneOwnjrJIT5pcbNykiR1jslJktQ5JidJUueYnCRJneOECHWekyOkxcfKSZLUOSYnSVLnOKynecUhPmlxsHKSJHWOyUmS1DkmJ0lS55icJEmd44QIzVtOjpAWLisnSVLnmJwkSZ3jsJ4WBIf4pIXFykmS1DkmJ0lS5zispwWnd4gPHOaT5iMrJ0lS55icJEmdY3KSJHWOyUmS1DlOiNCC53egpPnHykmS1DkmJ0lS5zisp0XFIT5pfhha5ZTkpCTXJFnTs+7oJJcluTjJZ5NsOdN92/WfSrK6vV2RZPWw4pckjc8wh/VOBpZPWncOsEdVPRj4PvCmDdiXqnpOVS2rqmXAZ4DT5yxaLTort9v79pukbhlacqqq84DrJq07u6rWt4tfB3ac6b69kgR4NnDa3EQrSeqScU6IeAlw1iz3fSzws6q6fKoGSY5IsirJqpWnnzrLh5EkjcNYJkQkeQuwHpht1jiUaaqmqloBrAA448K1NcvH0SLhRAmpW0aenJK8EDgQ2L+qNjhpJFkCPBN42FzHJknqhpEmpyTLgTcAj6uqm2d5mCcAl1XV2rmLTJIEsO7KnWe349ZzG8cwp5KfBpwP7J5kbZKXAscAdwPOaaeDf7htu32SM6fZd8IhOBFCQ+QsPmn8hlY5VdWhfVafOEXbq4EDptl3YtuLNjo4SVKnefkiSVLnmJwkSZ1jcpIkdY4XfpUG8PtP0nhYOUmSOsfkJEnqHIf1pBlyiE8aHSsnSVLnmJwkSZ3jsJ40C5MvbeQwnzS3rJwkSZ1jcpIkdY7DetIsbLfzlb+zfNHO295+/yHf+Nmow5EWHCsnSVLnWDlJMzS5WprKRQ/ftu96Kypp5qycJEmdY3KSJHWOw3rSADMdypuJ3uE+h/ikwaycJEmdY3KSJHWOw3rSGDjEJw1m5SRJ6hwrJ2nM/F6UdEdWTpKkzjE5SZI6x2E9qaMmD/c5zKfFxMpJktQ5JidJUuc4rCfNE343SouJlZMkqXNMTpKkznFYT5qHHOLTQmflJEnqHCsnaYB1V+58+/25/G2nuWQVpYVoaJVTkpOSXJNkTc+6o5NcluTiJJ9NsuUU+x6ZZE2SS5K8pmf9UUl+kmR1eztgWPFLksZnmMN6JwPLJ607B9ijqh4MfB940+SdkuwBvAzYB3gIcGCS+/U0eX9VLWtvZw4lcknSWA1tWK+qzkuydNK6s3sWvw4c3GfXBwBfr6qbAZJ8FTgIeM9wIpUWDq9wroVinBMiXgKc1Wf9GmDfJPdKsjlwALBTz/ZXt8OCJyXZahSBSpJGa2BySnL3JLv0Wf/gjXnQJG8B1gOnTt5WVZcC76YZAlwJXNS2BTgO2AVYBqwD3jvgMY5IsirJqpWn3+FhJEkj1G8ewiBTJqckzwYuAz7TTkzYu2fzyRsR4AuBA4HDqqr6tamqE6tqz6raF7gOuLxd/7OqurWqbgOOpzkv1VdVraiqvapqr+XPPGy24UoLwkUP3/b2mzQmJ3PHeQhTGlQ5vRl4WFUtA14MfCzJM9ttmU1kSZYDbwCeNnFOaYp2927/3Rl4JnBau7xdT7ODaIYAJUkdV1Xn0RQbMzJoQsSSqlrXHvSbSf4IOCPJjkDfiqdXktOA/YCtk6wF3kYzO+/3gHOSQDPx4RVJtgdOqKqJqeGfSXIv4BbgVVV1fbv+PUmWtY9/BfDymXZU2ljz4TtPM+H3ojQfDEpONybZpar+C6Cq1iXZD/gc8IfTHbiqDu2z+sQp2l5NM/FhYvmxU7R7/nSPK0kavSRHAEf0rFpRVStme7xByekNTBq+q6pftUNzd/h+kiRp8WoT0ayT0WSDktMpwD8meW9VrQdIsi3NDLndgXfMVRCSxsOfgldXDZoQ8TDgvsC3kzw+yZHAN4HzgYePIjhJ0sLQzkM4H9g9ydokLx3UfsrKqZ2E8Io2Kf0bcDXwiKpaO5cBS5IWvinmIUxpyuTUXpT13TRV0nKaCQtnJTmyqr68UVFK6iRn8qkrBp1zuhA4lmYq93rg7HYa97FJfryhWVCSpJkalJz2nTyEV1WrgUcledlww5IkLWaDzjlNeW6pqo4fTjjS/ND7hVyY31/KnYpDfBonf6ZdktQ5/ky7pGlZRWnUrJwkSZ1jcpIkdY7DetIcWChXLJ8Jfwpeo2DlJEnqHJOTJKlzHNaTNCdWbrf37feXr/vWGCPRQmDlJEnqHCsnaY4tpskRvXr7upK9p2xnVaWZsHKSJHWOyUmS1DkO60lD5BDfHS+S68QJzYSVkySpc0xOkqTOcVhP0tj0DvH1crhPVk6SpM4xOUmSOsdhPWlEnLnXmDx7rx+H+2TlJEnqHCsnSSM16DtQ07GiWjysnCRJnWNykiR1jsN60hgs1skRwzLVcB845DdfWTlJkjrH5CRJ6hyH9aQxW8xDfBszc2+mnOE3Pw21ckpyUpJrkqzpWffOJBcnWZ3k7CTb99nvj9rtE7f/TvKMSW0+mOSmYcYvSRqPYVdOJwPHAB/tWXd0Vf01QJI/B94KvKJ3p6r6CrCsbXNP4AfA2RPbk+wFbDnMwKVxmFw9LLZKapQGTaKYYHU1PkOtnKrqPOC6Setu7Fm8K1DTHOZg4KyquhkgySbA0cDr5zBUqXNMTFrMxjIhIsnfJLkKOIymchrkEOC0nuVXA1+oqnXTPMYRSVYlWbXy9FM3LmBJ0kilarrCZSMfIFkKnFFVe/TZ9ibgLlX1tin23Q64GNi+qm5pz099GtivqtYnuamqtpguhjMuXDvcTkpDsJgrp2FNjhiWcQ7/LXnKQZnL4832/fLAPXec0zjGPZX8E8CzBmx/NvDZqrqlXX4osCvwgyRXAJsn+cFwQ5QkjdrIk1OS+/UsPg24bEDzQ+kZ0quqL1bV71fV0qpaCtxcVbsOJ1JJ0rgMdbZektOA/YCtk6wF3gYckGR34Dbgx7Qz9doZeK+oqsPb5aXATsBXhxmj1FXO3Js/nPk394aanKrq0D6rT5yi7Srg8J7lK4Adpjn+tOebJEnzz7jPOUmSdAdevkiaJxbTZY5GcVmjUZvJ0N9k4xgKnP1j7jincVg5SZI6x8pJmocWUxW1mM2k2jpwBHGMg5WTJKlzTE6SpM5xWE+a5xb6EN9CnByh6Vk5SZI6x+QkSeock5MkqXNMTpKkznFChLSALPTJEVo8rJwkSZ1jcpIkdY7DetICtRCH+PzO0+Jh5SRJ6hyTkySpcxzWkxaBhTjEp4XNykmS1DkmJ0lS5zisJy0yDvFpPrBykiR1jpWTtIhN/q6QlZS6wspJktQ5JidJUuc4rCfpdlNdEsjhPo2alZMkqXNMTpKkznFYT9K0/G6URs3KSZLUOSYnSVLnOKwnaYM4xKdRsHKSJHWOlZOkWfN7URqWoVVOSU5Kck2SNX22vTZJJdl6in1vTbK6vX2hZ/2pSb6XZE17/E2HFb8kaXyGOax3MrB88sokOwFPBAZ9tPpNVS1rb0/rWX8qcH/gQcBmwOFzF64kqSuGNqxXVeclWdpn0/uB1wOfn8Uxz5y4n+SbwI6zjU/S8Ew13AcbN+Q36LhaWEY6ISLJ04CfVNVF0zS9S5JVSb6e5Bl9jrMp8Hxg5TDilCSN18iSU5LNgbcAb51B852rai/gucAHkuwyafuxwHlV9e8DHu+INsGtWnn6qbOOW5I0eqOcrbcLcF/goiTQDMldmGSfqvppb8Oqurr994dJzgUeCvwXQJK3AdsALx/0YFW1AlgBcMaFa2tOeyJp1hya00yMLDlV1XeAe08sJ7kC2Kuqru1tl2Qr4Oaq+m07m+/RwHvabYcDTwL2r6rbRhW7JGm0hjmV/DTgfGD3JGuTvHRA272SnNAuPgBYleQi4CvAu6rqu+22DwPbAue308xnMkQoSZpnhjlb79Bpti/tub+Kdlp4Vf0nzVTxfvv4pWFJWgS8fJEkqXNMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzjE5SZI6x+QkSeock5MkqXNMTpKkzhlLckpyUpJrkqzpWXdUkp8kWd3eDuiz312SfDPJRUkuSfL20UYuSRqFcVVOJwPL+6x/f1Uta29n9tn+W+DxVfUQYBmwPMkjhhinJGkMxpKcquo84LpZ7FdVdVO7uGl7q7mMTZI0fl075/TqJBe3w35b9WuQZJMkq4FrgHOq6hujDVGSNGxdSk7HAbvQDNetA97br1FV3VpVy4AdgX2S7NGvXZIjkqxKsmrl6acOK2ZJ0hAsGXcAE6rqZxP3kxwPnDFN+xuSnEtz7mpNn+0rgBUAZ1y41qE/SZpHOlM5JdmuZ/Eg+iScJNsk2bK9vxnwBOCy0UQoSRqVsVROSU4D9gO2TrIWeBuwX5JlNBMcrgBe3rbdHjihqg4AtgNOSbIJTWL9dFUNrLAkSfPPWJJTVR3aZ/WJU7S9GjigvX8x8NAhhiZJ6oDODOtJkjTB5CRJGroky5N8L8kPkrxxuvYmJ0nSULXzBD4EPBl4IHBokgcO2sfkJEkatn2AH1TVD6vqf4BPAk8ftIPJSZI0bDsAV/Usr23XTakzX8IdpgP33DHQXDWi/XLuvGY/umUh9GMh9AHsx1xY8pSDMpv9khwBHNGzakVPH/odc+DFERZb5XTE9E3mBfvRLQuhHwuhD2A/xqaqVlTVXj233uS6FtipZ3lH4OpBx1tsyUmSNHrfAu6X5L5J7gwcAnxh0A6LYlhPkjQ+VbU+yauBfwU2AU6qqksG7bPYktO8H4tu2Y9uWQj9WAh9APvRWe0PyPb7Edm+UuUFuyVJ3eI5J0lS5yyI5NT+cu41Se7wMxvt9qe3v7C7uv0Bwsf0bNs5ydlJLk3y3SRLRxV3nzgH9qOn3d5Jbk1y8KT1d0/ykyTHDDfSgbFN91rsl+SX7WuxOslbZ7rvKG1kP45MsibJJUleM7qo+8Y57XPa9mV1G+9Xe9ZvmeSfk1zW/n08cjRR941xutfjdT2vxZr27+OePds3SfLtJGP9FYMZ9OMeSf4lyUXt6/Hinm0rk9ww7j6MTFXN+xuwL7AnsGaK7Vvwv0OYDwYu69l2LvDEnnabd7UfbZtNgC/TjN0ePGnb3wOfAI7pah9ofirljNn2v+v9APag+S2yzWnO6f4bcL8O92NL4LvAzu3yvXu2nQIc3t6/M7BlV/sxqe1TgS9PWvd/2r+Nvv/3utIP4M3Au9v72wDXAXdul/dv+zbWPozqtiAqp6o6j+ZFnGr7TdW+usBdab/81V7baUlVndPT7uZhxzuV6frR+jPgM8A1vSuTPAzYFjh7ONHNzAz7MOf7zrWNiOUBwNer6uaqWg98lebHM8diBv14LnB6VV3Ztr8Gmiqc5o30xHb9/1TVDUMOd0ob+HocCpw2sZBkR+ApwAlDCG2DzKAfBdwtSWg+LF8HrG/3/RLwq6EH2RELIjnNRJKDklwGfBF4Sbt6N+CGJKe3Jf/R7QUKOynJDjRvdB+etP5OwHuB140jrll4ZDtscVaSPxx3MBuhXz/WAPsmuVeSzWl+i2ynqQ8xdrsBWyU5N8kFSV7Qrv8D4OfAR9q/jROS3HV8Yc5M+5wvp/kAN+EDwOuB28YS1IY5huYDztXAd4Ajq2o+xD3nFk1yqqrPVtX9gWcA72xXLwEeC7wW2JvmD/JFYwlwZj4AvKGqbp20/k+BM6vqqj77dM2FwH2q6iHAB4HPjTme2erbj6q6FHg3cA6wEriI9pNvRy0BHkZTWTwJ+Osku7Xr9wSOq6qHAr8Gpv2Zgw54KvAfVXUdQJIDgWuq6oLxhjVjTwJWA9sDy4Bj2ip20Vk0yWlCW1bvkmRrmktqfLuaK+Wup3mD2XOsAQ62F/DJJFcABwPHJnkG8Ejg1e36vwNekORdY4tygKq6sapuau+fCWzavhbzyqB+VNWJVbVnVe1LMyxz+RhDnc5aYGVV/bqqrgXOAx7Srl9bVd9o2/0z3f7bmHAIPUN6wKOBp7V/G58EHp/k4+MIbIZeTDPMWlX1A+BHwP3HHNNYLIrklGTXdgyXJHvSnNz9Bc0lNbZKsk3b9PE0J4c7qaruW1VLq2opzZvFn1bV56rqsKrauV3/WuCjVdXJT7lJfr/ntdiH5v/gL8Yb1YYb1I8k927/3Rl4Jr/7Ztk1nwcem2RJOyT2cODSqvopcFWS3dt2+9Phvw1oZroBj6PpEwBV9aaq2rH92ziEZqLE88YU4kxcSfNck2RbYHfgh2ONaEwWxBUikpxGM3tq6yRrgbcBmwJU1YeBZ9FUE7cAv6s04TwAAAGqSURBVAGe006QuDXJa4EvtW80FwDHj6ELwIz60Xkz6MPBwCuTrKd5LQ6ZmKzSb9+qOnHknZgiFmbYD+AzSe4F3AK8qqquH3X8E6brR1VdmmQlcDHNOZkTqmpimvOfAaemuRbaD2k+1Y/FDP82DgLOrqpfjyXIGZhBP94JnJzkOzRX8n5DW9GS5N9pqqgt2n1fWlX/OvpejIZXiJAkdc6iGNaTJM0vJidJUueYnCRJnWNykiR1jslJktQ5JidplpLslORHE1e/TrJVu3yfRXcFaWmOmZykWWovF3UcMHE1jncBK6rqx8DRwPPHFZs035mcpI3zfuARaX636TE0F+BddFeQlubagrhChDQuVXVLktfRXOT1j6vqf8Ydk7QQWDlJG+/JwDqaHxqUNAdMTtJGSLIMeCLwCOAvkmw35pCkBcHkJM1Se7Hg44DXtL8kezTNT5ZI2kgmJ2n2XgZcWVXntMvHAvdP8rj2CtL/BOyfZG2SJ40tSmke8qrkkqTOsXKSJHWOyUmS1DkmJ0lS55icJEmdY3KSJHWOyUmS1DkmJ0lS55icJEmd8/8BYST1VBDrAaIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "qda_2d = QDA()\n",
    "qda_2d.fit(X_2d, y)\n",
    "graph_boundaries(X_2d, qda_2d, 'QDA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive Bayes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we implement a Naive Bayes model below. This model allows us to assign each variable in our dataset a distribution, though by default they are all assumed to be Normal. Since each variable has its own distribution, estimating the model's parameters is more involved. For each variable and each class, we estimate the parameters separately through the `_estimate_class_parameters`. The structure below allows for Normal, Bernoulli, and Poisson distributions, though any distribution could be implemented.\n",
    "\n",
    "Again, we make predictions by calculating $p(Y_n = k|\\bx_n)$ for $k = 1, \\dots, K$ through Bayes' rule and predicting the class with the highest posterior probability. Since each variable can have its own distribution, this problem is also more involved. The `_get_class_probability` method calculates the probability density of a test observation's input variables. By the conditional independence assumption, this is just the product of the individual densities. \n",
    "\n",
    "Naive Bayes performs worse than LDA or QDA on the training data, suggesting the conditional independence assumption might be inappropriate for this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NaiveBayes:\n",
    "    \n",
    "    ######## Fit Model ########\n",
    "\n",
    "    def _estimate_class_parameters(self, X_k):\n",
    "        \n",
    "        class_parameters = []\n",
    "        \n",
    "        for d in range(self.D):\n",
    "            X_kd = X_k[:,d] # only the dth column and the kth class\n",
    "            \n",
    "            if self.distributions[d] == 'normal':\n",
    "                mu = np.mean(X_kd)\n",
    "                sigma2 = np.var(X_kd)\n",
    "                class_parameters.append([mu, sigma2])\n",
    "            \n",
    "            if self.distributions[d] == 'bernoulli':\n",
    "                p = np.mean(X_kd)\n",
    "                class_parameters.append(p)\n",
    "                \n",
    "            if self.distributions[d] == 'poisson':\n",
    "                lam = np.mean(X_kd)\n",
    "                class_parameters.append(p)\n",
    "                \n",
    "        return class_parameters\n",
    "    \n",
    "    def fit(self, X, y, distributions = None):\n",
    "        \n",
    "        ## Record info\n",
    "        self.N, self.D = X.shape\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        if distributions is None:\n",
    "            distributions = ['normal' for i in range(len(y))]\n",
    "        self.distributions = distributions\n",
    "        \n",
    "        \n",
    "        ## Get prior probabilities \n",
    "        self.unique_y, unique_y_counts = np.unique(self.y, return_counts = True) # returns unique y and counts\n",
    "        self.pi_ks = unique_y_counts/self.N\n",
    "        \n",
    "        \n",
    "        ## Estimate parameters\n",
    "        self.parameters = []\n",
    "        for i, k in enumerate(self.unique_y):\n",
    "            X_k = self.X[self.y == k]\n",
    "            self.parameters.append(self._estimate_class_parameters(X_k))\n",
    "    \n",
    "    \n",
    "    ######## Make Classifications ########\n",
    "            \n",
    "    def _get_class_probability(self, x_n, j):\n",
    "        \n",
    "        class_parameters = self.parameters[j] # j is index of kth class\n",
    "        class_probability = 1 \n",
    "        \n",
    "        for d in range(self.D):\n",
    "            x_nd = x_n[d] # just the dth variable in observation x_n\n",
    "            \n",
    "            if self.distributions[d] == 'normal':\n",
    "                mu, sigma2 = class_parameters[d]\n",
    "                class_probability *= sigma2**(-1/2)*np.exp(-(x_nd - mu)**2/sigma2)\n",
    "            \n",
    "            if self.distributions[d] == 'bernoulli':\n",
    "                p = class_parameters[d]\n",
    "                class_probability *= (p**x_nd)*(1-p)**(1-x_nd)\n",
    "                \n",
    "            if self.distributions[d] == 'poisson':\n",
    "                lam = class_parameters[d]\n",
    "                class_probability *= np.exp(-lam)*lam**x_nd\n",
    "                \n",
    "        return class_probability \n",
    "            \n",
    "    def classify(self, X_test):\n",
    "        \n",
    "        y_n = np.empty(len(X_test))\n",
    "        for i, x_n in enumerate(X_test): # loop through test observations\n",
    "            \n",
    "            x_n = x_n.reshape(-1, 1)\n",
    "            p_ks = np.empty(len(self.unique_y))\n",
    "        \n",
    "            for j, k in enumerate(self.unique_y): # loop through classes\n",
    "                    \n",
    "                p_x_given_y = self._get_class_probability(x_n, j)\n",
    "                p_y_given_x = self.pi_ks[j]*p_x_given_y # bayes' rule\n",
    "\n",
    "                p_ks[j] = p_y_given_x\n",
    "            \n",
    "            y_n[i] = self.unique_y[np.argmax(p_ks)]\n",
    "        \n",
    "        return y_n\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9775280898876404"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nb = NaiveBayes()\n",
    "nb.fit(X, y)\n",
    "yhat = nb.classify(X)\n",
    "np.mean(yhat == y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAFNCAYAAACkMKB8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debhkVXnv8e9PGoMzKGhAaFEQHBAaRcQ4EZDYoqgYY8BZQTQJBnNFcYo4XL0qMeoV0TSDOCDqFRyC2EJUROMALTbYCCpBhBYEFRAQJ+S9f+x1sDjUGbr71Knd3d/P89Rz9rD23u+qqlNvrbXX3pWqQpKkPrnduAOQJGkyk5MkqXdMTpKk3jE5SZJ6x+QkSeodk5MkqXdMTuuQJK9Ncsy441jXJHljko+16YVJbkiywWrsZ95fnySVZNv5POZ8S7J7kpXjjkNzy+TUI0kuSXJlkjsNLDswyRmz2b6q3lZVB44grjOS/K59KP86yZlJHjLXx1ld7cPp5hbf9Ul+mOSFozhWVV1aVXeuqj/NIqZbfWCO6vWZD0lelmRFktsPLHt5ku8lWdDml7Tn/uYkLxhbsAOSbJnkmiSPHli2VVv2iDb/b0l+3N47FyZ53vgi1gSTU/8sAA4ZdxBDHFxVdwbuAZwBfHS84dzG5S2+uwKHAUcnedDkQhMfpFpl7weuBV4HkOR+wJuAA6rqplbmXOAfgXPGEuEQVbWS7v1wTJKN2uL/AD5UVd9p878B9gHuBjwfeG+Sv5r3YHUrJqf+OQI4NMnGw1YmeW+Sy5Jcl+S7SR4zsG6w+2lpkoMnbXtukqe36QckOT3J1e3b7jNnE1z7IPoEcMsHf5Jdk3wrybVJrkhy5MQ37CTvT/KuSXH8Z5KXt+ktkpyU5BdJfpLknyftd1mr65VJ/n0W8VVVfRa4BnhQkq1b19YBSS4FvtL2vVuSb7aYz02y+8Bx75vka+2b9OnApgPrJvY30Vq4e5IPJbm8fRv/bGv5fhHYorXmbmj1vOX1ads+Jcn5LYYzkjxwYN0lSQ5Ncl5rrX5y4sM1yaZJTmnbXZ3k60mm+1/eO8nFSX6Z5Igkt0vyF23bW1rASe6Z5LdJNhvyvN4MHAD8S5IdgaOBo6rqnIEy76+qLwO/m+l1SvKk1uq6rr2f3zjkOX5+kktb3K8bWH+HJMe35/sHwMNnONzRwBXA4UmeD2wPvH4g7sOr6sKqurklrK8Dj5ypDhqxqvLRkwdwCfB44GTgf7dlBwJnDJR5Dl3rZQHwCuDnwEZt3RuBj7Xp5wH/PbDdg+i++f4FcCfgMuCFbT8PBX4JPHiKuM4ADmzTtwfeCpw5sP5hwG5tX1sDFwAvb+t2BS4HbtfmNwVuBO5F9+Xou8Ab2n7vB1wMPKGV/Rbw3DZ9Z2C3KeLbHVjZpm8H7Av8ke5DaGuggI+0et8BuDfwK2DvVn6vNr/ZwHH/vT1XjwWuH3heJ/a3oM1/AfgksAmwIfC4yTENxDn4+mxH9419r7bdq4CLgNsPvBfOArYA7t6e05e2df8H+GDbbkPgMUCmeG4K+Grbx0LgRwOv5VHAOwbKHgL85wzv0de098oPae+7IWW+Abxghv3sDjykPf87AlcCT5v0HB/dXq+dgN8DD2zr306XQO4ObAWsmPxcDzneNsCv6b607DFNuTvQJbLF4/48WN8fYw/Ax8CL8efktEP7R9qMSclpyDbXADu16cEPv7u0D7/7tPm3Ase16b8Hvj5pP/8BHD7FMc6gSyjXAn9ose05TUwvBz4zMH8BsFebPhg4tU0/Arh00ravoetyATiTruto0xmet92Bm1t8VwPLgf3auokPuvsNlD8M+OikfXyJrktnIXATcKeBdR9nSHICNm/H3WSKmKZLTv8KfGpg3e2AnwG7D7wXnjOw/p3AB9v0m4HPAdvO4j1Vgx+0dN1uXx54/i/jz18clgHPnGF/j277fOs0ZWZMTkO2eQ/w7knP8ZYD688aeE0vnlSngyY/10P2v6C9D39K+2IxRbkPA0uZItn7mL+H3Xo9VFUrgFOAV09el+QVSS5oXT3X0vWTbzpkH9fTfavfry3aDzihTd8HeETrFrq27efZwF9OE9Y/V9XGwEbAk4FPt+4dkmzXupl+nuQ64G2TYvowXYuP9nfifNV96Lq+BuN4LV2rCrpupO2AC5OcneTJ08R3eVVtXFV3r6pFVfWJSesvG5i+D/B3k477aLpkswVwTVX9ZqD8T6c45lbA1VV1zTRxTWWLwf1W1212GV2rbsLPB6ZvpGs9Qtf1exFwWuuuu837ZJLBuv+0HZvqurB+AzwuyQOAbYHPT7WT1lX7H8D7gIPTnXdaLUkekeSrrTv318BLue37eKr6bzGkTjN5NV3r+Crg0CliOoLui+Ezq2UqjY8nh/vrcLoTy7ecr0l3fukwYE/g/Kq6Ock1QKbYx4l0/exn0nVXfLUtvwz4WlXttapBtQ/Rrye5CPgb4DzgA8D3gP2r6vp2PukZA5t9DFiRZCfggcBnB+L4SVXdf4pj/RjYv51PeTpdQrzHpMQx69AHpi+jazm9eHKhJPcBNklyp4HjLJy0/eB+7p5k46q6dprjDXM5XbfWxHFDl+x+NsN2E188XgG8IsmDga8mObu68z3DbAWc36YXtmNPmPji8HPg01U13fmif6X7cD8E+C1dolrl91DzceBI4IlV9bsk72HIl6wpXMFt6zSldANjXknXUrw98I0kJ7X310SZNwFPpOuWvW6VaqKRsOXUU1V1Ed25jH8eWHwXui6nXwALkryBbnTaVE6layW8GfhkSyzQtcq2S/LcJBu2x8MHT8hPJ8kj6c5hTXw43AW4DrihfQP/h0l1WQmcTddiOqmqfttWnQVcl+SwdpJ7gyQ7JHl4O85zkmzW4p748J92CPcsfQzYJ8kT2jE3Sjf0e8uq+ild99abktw+3RDkfYbtpKquoBv4cFSSTdrz+Ni2+krgHknuNkUMnwKelGTPJBvSJZvfA9+cKfgkT06ybUto19E9J9M9L69s8W1Fl1g+ObDuo3Tn6J5Dd15uqmPuRPdefHFrVbwR2DoDQ/bb87UR3ZelDdvzOtVnzF3oWp2/S7Ir8Kzpa30rnwJe0+q0JfCyaeK+HXAs8M7qBj2cB/xfYEl7/kjymnb8varqV6sQh0bI5NRvb6Y7iT/hS3Qfhj+i68r4Hbfu3riVqvo93eCKx9N9U51Yfj1dq2c/um/RPwfeQTcAYCpHpo08o/tAe31VfbGtO5Tun/t6upPYnxyy/YfpWgq3DEGv7lqhfYBFwE/oTrQfQ9dVCbAYOL8d87105xxmHAk2k6q6DHgqXRfiL+iew1fy5/+HZ9F9y76argU75Yc28Fy6wRcX0rUqXt6OcSFdy/Xi1nW4xaQYfkiXEN5HV+99gH2q6g+zqML9gf8CbqAbvHFUVZ0xTfnP0Q08WU7X1XvsQBwr6VroRTfI4DbSXXB8LN15povadr8FXgwckWSiG/Y0uhbVXwFL2vRjb7tHoDv39eYk19MNiPnUtDW+tTfRvf9/0o453WUNhwB3pDtnN+EtdF3YE9ecvY2u9fXj/Hl05WtXIR6NQOxa1XxoLYqPAVsPtODUA0mOoztn9/oZC0vzxHNOGrnWbXUIcIyJqV+SbE13Pm/n8UYi3Zrdehqpdh7rWrqRcO8ZczgakOQtdNcIHVFVPxl3PFp3pbtl1FfbSOPzk8x4Fxy79SRJI5Vkc2DzqjonyV3ozoE+rap+MNU2Y2k5JTkuyVVJVgwsu3u62+n8uP3dZIptFyY5rWXgH7RuCUlST1XVFdVuddUGZF3Ara/pu41xdesdTzcSa9Cr6a5cvz/wZYZcgNp8hK4b4oF0t8a5alRBSpLmVmtQ7Ax8Z7pyYxkQUVVnDmnxPJXuli/QDTs+g+6C01u0i+kWVNXpbT83zOZ4p5yz0r7L9cjiK84edwjSvFnwpH2nugh/tXz3l99crc/LXTZ71EvobiU1YUlVLRksk+TOwEl0996c9mLnPg2IuFe7qHHi4sZ7DimzHXBtkpPT3dH4iEzxo29JDkp3R+tlS08+YVgRSdIcqaolVbXLwGNyYtqQLjGdUFUnz7S/tW0o+QK6OzDvDFxKd7HnCxi4qHBCe2KWgC2n9YGtJam/2t04jgUuqKoZf/oG+tVyurKN6JgY2THsXNJK4HtVdXF1vyv0Wbqfe5Ak9dej6O6mskeS5e2x93Qb9Knl9Hm6nyx4e/v7uSFlzqa7KedmVfULYA+6+6BJknqqqr7B1DeoHmpcQ8lPpLsn2PZJViY5gC4p7ZXkx3R3On57K7tLkmPglnuxHQp8Ocn36Sp79DjqIEkanXGN1tt/ilV7Dim7jD/foJE2Um/HEYUmSeqBPp1zkiQJMDlJknrI5CRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTe6dO99aRV4p3IpXWXLSdJUu+YnCRJvWNykiT1jslJktQ7JidJUu+YnCRJvWNykiT1jslJktQ7JidJUu+YnCRJvWNykiT1jslJktQ7JidJUu+YnCRJvWNykiT1jslJktQ7JidJUu+YnCRJvWNykiT1jslJktQ7JidJUu+YnCRJvWNykiT1jslJktQ7I01OSY5LclWSFQPL/i7J+UluTrLLFNttlOSsJOe2sm8aWLdnknOSLE/yjSTbjrIOkqT5N+qW0/HA4knLVgBPB86cZrvfA3tU1U7AImBxkt3aug8Az66qRcDHgdfPacSSpLFbMMqdV9WZSbaetOwCgCTTbVfADW12w/aoidXAXdv03YDL5yxgSVIvjDQ5rYkkGwDfBbYF3l9V32mrDgROTfJb4Dpgtyl2IUlaS/V2QERV/al13W0J7Jpkh7bqX4C9q2pL4EPAvw/bPslBSZYlWbb05BPmJ2hJ0pzobctpQlVdm+QMuvNOVwI7DbSiPgksnWK7JcASgFPOWVnDymjttnTzh98yvfiKs8cYiaS51suWU5LNkmzcpu8APB64ELgGuFuS7VrRvYALxhOlJGlURtpySnIisDuwaZKVwOHA1cD7gM2ALyRZXlVPSLIFcExV7Q1sDny4nXe6HfCpqjql7fPFwElJbqZLVi8aZR0kSfNv1KP19p9i1WeGlL0c2LtNnwfsPMU+PzNse0nSuqOX3XqSpPWbyUmS1DsmJ0lS75icJEm9Y3KSJPWOyUmS1DsmJ0lS75icJEm9Y3KSJPVO72/8Ks2GN4GV1i22nCRJvWNykiT1jslJktQ7JidJUu+YnCRJvWNykiT1jslJktQ7Xuekdc7gNU/gdU/S2siWkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpdxytJ0lrmcERqU8eYxyjZMtJktQ7JidJUu/Yrad1nj9EKK19bDlJknrH5CRJ6h2TkySpd0xOkqTecUCE1isOjtDaavLd9td1tpwkSb1jcpIk9Y7depLUU+tbV96gkbWckhyX5KokKwaWHZHkwiTnJflMko1nu21b/skky9vjkiTLRxW/JGl8RtmtdzyweNKy04EdqmpH4EfAa1ZhW6rq76tqUVUtAk4CTp6zaLXeWbr5w295SOqXkSWnqjoTuHrSstOq6qY2+21gy9luOyhJgGcCJ85NtJKkPhnngIgXAV9czW0fA1xZVT+eqkCSg5IsS7Js6cknrOZhJEnjMJYBEUleB9wErG7W2J8ZWk1VtQRYAnDKOStrNY+j9YTXP6kP7GL+s3lPTkmeT/f7WHtW1SonjSQLgKcDD5vr2CRJ/TCvySnJYuAw4HFVdeNq7ubxwIVVtXLuIpMkAVxx6cLV23DTuY1jlEPJTwS+BWyfZGWSA4AjgbsAp7fh4B9sZbdIcuoM207YDwdCSNI6bWQtp6raf8jiY6coezmw9wzbTqx7wRoHJ0nqNW9fJEnqHW9fJE3iyD3NJ0foDWfLSZLUO7acpGnYitIo2FqamS0nSVLvmJwkSb1jt540S3bxaU3YlbdqbDlJknrH5CRJ6h279aTVMLmLxm4+DWNX3uqz5SRJ6h2TkySpd+zWk+aAI/k0wa68uWHLSZLUO7acpDlmK2r9Y2tp7tlykiT1jslJktQ7dutJI2QX37rLrrzRsuUkSeodk5MkqXfs1pPmiV18aze78eaXLSdJUu/YcpLGwFZUf9lC6gdbTpKk3jE5SZJ6x249acym6kayu2/+2JXXP7acJEm9Y3KSJPWO3XpSj2y+8NJbps9deK+RH2+n71w58mP0id13aw9bTpKk3jE5SZJ6x249acwGu/Lm27mPWPWuw7WhK9Duu7WfLSdJUu/YcpLGYJytpTU1m9bWXLaubAWtn0bWckpyXJKrkqwYWHZEkguTnJfkM0k2nmLbQ5KsSHJ+kpcPLH9jkp8lWd4ee48qfknS+IyyW+94YPGkZacDO1TVjsCPgNdM3ijJDsCLgV2BnYAnJ7n/QJF3V9Wi9jh1JJFLksZqZN16VXVmkq0nLTttYPbbwDOGbPpA4NtVdSNAkq8B+wLvHE2kkubadF1/V1y6cB4j0dpqnAMiXgR8ccjyFcBjk9wjyR2BvYGtBtYf3LoFj0uyyXwEKkmaX9MmpyR3TbLNkOU7rslBk7wOuAk4YfK6qroAeAddF+BS4NxWFuADwDbAIuAK4F3THOOgJMuSLFt68m0OI0maR8PGIUxnym69JM8E3gNclWRD4AVVNXGb5OOBh65mgM8HngzsWVU1rExVHQsc28q/DVjZll85sJ+jgVOmOk5VLQGWAJxyzsqhx5Hm09o8Qm8uDT4PdvGtV44HjgQ+MpvC07WcXgs8rKoWAS8EPprk6W1dVieyJIuBw4CnTJxTmqLcPdvfhcDTgRPb/OYDxfal6wKUJPVcVZ0JXD3b8tMNiFhQVVe0nZ6V5K+BU5JsCczYEklyIrA7sGmSlcDhdKPz/gI4PQl0Ax9emmQL4JiqmhgaflKSewB/BP6pqq5py9+ZZFE7/iXAS2ZbUWkcbC1Jq2e65HRdkm2q6n8AquqKJLsDnwUePNOOq2r/IYuPnaLs5XQDHybmHzNFuefOdFxJ0vxLchBw0MCiJe30ymqZLjkdxqTuu6q6vnXN3eb6JEnS+mvwPP9cmO6c04eBv01ySwJLci/gQ8A+cxWApPXX5gsvveUhDZouOT0MuC/wvSR7JDkEOAv4FvCI+QhOkrRuaOMQvgVsn2RlkgOmKz9lt14bhPDSlpT+C7gc2K2qVs5lwJKkdd8U4xCmNN11ThvTXQz7CLp75O0NfDHJIVX1lTWKUlqH2UUlrbnpBkScAxxFN5T7JuC0Noz7qCQ/XdUsKEnSbE2XnB47uQuvqpYDf5XkxaMNS5K0PptyQMR055aq6ujRhCNJkj/TLknqIX+mXVIveENYDbLlJEnqHZOTJKl3TE6SpN4xOUmSesfkJEnqHUfrSXPAWxZJc8uWkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrH65wk9Y53KJctJ0lS75icJEm9Y3KSJPWOyUmS1DsmJ0lS75icJEm9Y3KSJPWOyUmS1DsmJ0lS75icJEm94+2LJPWOtyzSSFtOSY5LclWSFQPL3pLkvCTLk5yWZIsh2/11Wz/x+F2Sp00q874kN4wyfknSeIy6W+94YPGkZUdU1Y5VtQg4BXjD5I2q6qtVtaiV2QO4EThtYn2SXYCNRxa1JGmsRpqcqupM4OpJy64bmL0TUDPs5hnAF6vqRoAkGwBHAK+aw1AlST0ylgERSd6a5DLg2QxpOU2yH3DiwPzBwOer6ooZjnFQkmVJli09+YQ1C1iSNK/Gkpyq6nVVtRVwAl2yGSrJ5sBDgC+1+S2AvwPeN4tjLKmqXapql8VPf/bcBC5JmhfjHkr+ceBvp1n/TOAzVfXHNr8zsC1wUZJLgDsmuWi0IUqS5tu8J6ck9x+YfQpw4TTF92egS6+qvlBVf1lVW1fV1sCNVbXtaCKVJI3LSK9zSnIisDuwaZKVwOHA3km2B24Gfgq8tJXdBXhpVR3Y5rcGtgK+NsoYpbkweF3O4E+MS1o9I01OVbX/kMXHTlF2GXDgwPwlwL1n2P+d1yQ+SVI/jfuckyRJt+HtiyRJt1h8xdmrueWWcxqHLSdJUu/YcpLUC97sVYNsOUmSesfkJEnqHbv1pDnmNU/SmrPlJEnqHZOTJKl37NaTNDaO0NNUbDlJknrHlpM0Qg6OkFaPLSdJUu+YnCRJvWO3nqR55SAIzYYtJ0lS75icJEm9Y7eeNE/W15F7duNpddhykiT1jslJktQ7dutJY7Cud/HZlac1ZctJktQ7tpykMVtXWlG2ljSXbDlJknrH5CRJ6h279aQeWdu6+OzK06jYcpIk9Y7JSZLUO3brST01XZfZfHf52X2n+WbLSZLUOyYnSVLv2K0nrYXsZtO6zpaTJKl3TE6SpN4ZWXJKclySq5KsGLLu0CSVZNMptv1TkuXt8fmB5Sck+WGSFW3/G44qfknS+Iyy5XQ8sHjywiRbAXsB042F/W1VLWqPpwwsPwF4APAQ4A7AgXMXriSpL0aWnKrqTODqIaveDbwKqNXY56nVAGcBW65ZlJKkPprXc05JngL8rKrOnaHoRkmWJfl2kqcN2c+GwHOBpaOIU5I0XvOWnJLcEXgd8IZZFF9YVbsAzwLek2SbSeuPAs6sqq9Pc7yDWoJbtvTkE1Y7bknS/JvP65y2Ae4LnJsEui65c5LsWlU/HyxYVZe3vxcnOQPYGfgfgCSHA5sBL5nuYFW1BFgCcMo5K1e5C1GSND7zlpyq6vvAPSfmk1wC7FJVvxwsl2QT4Maq+n0bzfco4J1t3YHAE4A9q+rm+YpdkjS/RjmU/ETgW8D2SVYmOWCasrskOabNPhBYluRc4KvA26vqB23dB4F7Ad9qw8xn00UoSVrLjKzlVFX7z7B+64HpZbRh4VX1Tbqh4sO28XZLkrQe8A4RkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknpnLMkpyXFJrkqyYmDZG5P8LMny9th7yHYbJTkryblJzk/ypvmNXJI0H8bVcjoeWDxk+buralF7nDpk/e+BPapqJ2ARsDjJbiOMU5I0BmNJTlV1JnD1amxXVXVDm92wPWouY5MkjV/fzjkdnOS81u23ybACSTZIshy4Cji9qr4zvyFKkkatT8npA8A2dN11VwDvGlaoqv5UVYuALYFdk+wwrFySg5IsS7Js6cknjCpmSdIILBh3ABOq6sqJ6SRHA6fMUP7aJGfQnbtaMWT9EmAJwCnnrLTrT5LWIr1pOSXZfGB2X4YknCSbJdm4Td8BeDxw4fxEKEmaL2NpOSU5Edgd2DTJSuBwYPcki+gGOFwCvKSV3QI4pqr2BjYHPpxkA7rE+qmqmraFJUla+4wlOVXV/kMWHztF2cuBvdv0ecDOIwxNktQDvenWkyRpgslJkjRySRYn+WGSi5K8eqbyJidJ0ki1cQLvB54IPAjYP8mDptvG5CRJGrVdgYuq6uKq+gPwCeCp021gcpIkjdq9gcsG5le2ZVPqzUW4o/Tkh24Z6O4a0S7OXatZj35ZF+qxLtQBrMdcWPCkfbM62yU5CDhoYNGSgToM2+e0N0dY31pOB81cZK1gPfplXajHulAHsB5jU1VLqmqXgcdgcl0JbDUwvyVw+XT7W9+SkyRp/p0N3D/JfZPcHtgP+Px0G6wX3XqSpPGpqpuSHAx8CdgAOK6qzp9um/UtOa31fdGN9eiXdaEe60IdwHr0VvsB2WE/IjtUqrxhtySpXzznJEnqnXUiObVfzr0qyW1+ZqOtf2r7hd3l7QcIHz2wbmGS05JckOQHSbaer7iHxDltPQbKPTzJn5I8Y9Lyuyb5WZIjRxvptLHN9FrsnuTX7bVYnuQNs912Pq1hPQ5JsiLJ+UlePn9RD41zxue01WV5i/drA8s3TvLpJBe2/49Hzk/UQ2Oc6fV45cBrsaL9f9x9YP0GSb6XZKy/YjCLetwtyX8mObe9Hi8cWLc0ybXjrsO8qaq1/gE8FngosGKK9Xfmz12YOwIXDqw7A9hroNwd+1qPVmYD4Ct0fbfPmLTuvcDHgSP7Wge6n0o5ZXXr3/d6ADvQ/RbZHenO6f4XcP8e12Nj4AfAwjZ/z4F1HwYObNO3Bzbuaz0mld0H+MqkZf+r/W8Mfe/1pR7Aa4F3tOnNgKuB27f5PVvdxlqH+XqsEy2nqjqT7kWcav0N1V5d4E60i7/avZ0WVNXpA+VuHHW8U5mpHs3LgJOAqwYXJnkYcC/gtNFENzuzrMOcbzvX1iCWBwLfrqobq+om4Gt0P545FrOox7OAk6vq0lb+Kuha4XQfpMe25X+oqmtHHO6UVvH12B84cWImyZbAk4BjRhDaKplFPQq4S5LQfVm+Gripbftl4PqRB9kT60Rymo0k+ya5EPgC8KK2eDvg2iQntyb/Ee0Ghb2U5N50H3QfnLT8dsC7gFeOI67V8MjWbfHFJA8edzBrYFg9VgCPTXKPJHek+y2yrabexdhtB2yS5Iwk303yvLb8fsAvgA+1/41jktxpfGHOTnvOF9N9gZvwHuBVwM1jCWrVHEn3Bedy4PvAIVW1NsQ959ab5FRVn6mqBwBPA97SFi8AHgMcCjyc7h/yBWMJcHbeAxxWVX+atPwfgVOr6rIh2/TNOcB9qmon4H3AZ8ccz+oaWo+qugB4B3A6sBQ4l/bNt6cWAA+ja1k8AfjXJNu15Q8FPlBVOwO/AWb8mYMe2Af476q6GiDJk4Grquq74w1r1p4ALAe2ABYBR7ZW7HpnvUlOE1qzepskm9LdUuN71d0p9ya6D5iHjjXA6e0CfCLJJcAzgKOSPA14JHBwW/5vwPOSvH1sUU6jqq6rqhva9KnAhu21WKtMV4+qOraqHlpVj6XrlvnxGEOdyUpgaVX9pqp+CeFVYg8AAAKfSURBVJwJ7NSWr6yq77Ryn6bf/xsT9mOgSw94FPCU9r/xCWCPJB8bR2Cz9EK6btaqqouAnwAPGHNMY7FeJKck27Y+XJI8lO7k7q/obqmxSZLNWtE96E4O91JV3beqtq6qrek+LP6xqj5bVc+uqoVt+aHAR6qql99yk/zlwGuxK9178FfjjWrVTVePJPdsfxcCT+fWH5Z98zngMUkWtC6xRwAXVNXPgcuSbN/K7UmP/zegG+kGPI6uTgBU1Wuqasv2v7Ef3UCJ54wpxNm4lO65Jsm9gO2Bi8ca0ZisE3eISHIi3eipTZOsBA4HNgSoqg8Cf0vXmvgj8Fvg79sAiT8lORT4cvug+S5w9BiqAMyqHr03izo8A/iHJDfRvRb7TQxWGbZtVR0775WYIhZmWQ/gpCT3AP4I/FNVXTPf8U+YqR5VdUGSpcB5dOdkjqmqiWHOLwNOSHcvtIvpvtWPxSz/N/YFTquq34wlyFmYRT3eAhyf5Pt0d/I+rLVoSfJ1ulbUndu2B1TVl+a/FvPDO0RIknpnvejWkyStXUxOkqTeMTlJknrH5CRJ6h2TkySpd0xO0mpKslWSn0zc/TrJJm3+PuvdHaSlOWZyklZTu13UB4CJu3G8HVhSVT8FjgCeO67YpLWdyUlaM+8Gdkv3u02PprsB73p3B2lprq0Td4iQxqWq/pjklXQ3ef2bqvrDuGOS1gW2nKQ190TgCrofGpQ0B0xO0hpIsgjYC9gN+Jckm485JGmdYHKSVlO7WfAHgJe3X5I9gu4nSyStIZOTtPpeDFxaVae3+aOAByR5XLuD9P8D9kyyMskTxhaltBbyruSSpN6x5SRJ6h2TkySpd0xOkqTeMTlJknrH5CRJ6h2TkySpd0xOkqTeMTlJknrn/wPoAkw7hIOGKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_2d = NaiveBayes()\n",
    "nb_2d.fit(X_2d, y)\n",
    "graph_boundaries(X_2d, nb_2d, 'Naive Bayes')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
